Projective studies of spin nematics in a quantum frustrated ferromagnet 
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We study the ground state properties of the spin-i frustrated ferromagnetic J1-J2 Heisenberg 
model on the square lattice, employing projected BCS wavefunctions with spin-triplet pairings of 
the spinon fields as trial wavefunctions. Based on the variational Monte Carlo analysis, we argue 
that, in the competing coupling regime, a certain type of the projected BCS wavefunction, dubbed 
the projected Z2 planar state, achieves the best optimal energy among the other competing states 
such as the ferromagnetic state and coUinear antiferromagnetic state. Like in quantum spin liquids, 
the projected Z2 planar state preserves the translational symmetry of the square lattice. However, it 
is also accompanied by a d-wave ordering of the quadrupole moments, breaking the spin rotational 
symmetry. The state thus describes a quantum spin analogue of the nematic liquid crystals. The 
calculated static correlation functions also reveal that the projected Z2 planar state has a strong 
coUinear antiferromagnetic fluctuation. 

PACS numbers: 



I. INTRODUCTION 

Deciphering unidentified states of quantum zero-point 
motion - quantum vacuum - is one of the central re- 
search field in condensed matter physics li In the realm 
of quantum magnetism, investigations of spin-rotational 
symmetric quantum spin liquid (QSL) - a quantum mag- 
net which remains totally disordered all the way down 
to the zero temperature - belong to this research cate- 
gory. Indeed, after seminal proposal of spin-singlet res- 
onating valence bond (RVB) wavefunctions by Anderson 
and his co-workers?^ tremendous research efforts are de- 
voted to establishing a true realization of QSL or RVB- 
type ground state in spatial dimension greater than onci^ 
Among others, 'frustrated' Mott insulating magnets are 
regarded as promising candidate materials for this in- 
vestigation)^ where competing magnetic interactions be- 
tween localized spins often make it hard for a system to 
fall into a simple classical spin ordering. 

This work reports a variational study of spin-i quan- 
tum frustrated ferromagnetic model. One motivation of 
this research is a couple of experimental works on two- 
dimensional films of solid ■^Hc;^— where the interactions 
between S — ^ nuclear spins of ^He atoms are highly frus- 
trated but predominantly ferromagnetic ii The specific 
heat measurement^ and magnetic susceptibility^ at the 
ultra-low temperature regime conclude that the ground 
state of this frustrated quantum spin system is a QSL-like 
state with either a gapless spin excitation or an extremely 
small spin gap. These experiments suggest a possibility 
of an exotic quantum phase in quantum frustrated ferro- 
magnets. 

The other incentive of this work stems from recent 
theoretical studies of frustrated magnets with ferromag- 



netic nearest-neighbor interactions on the square lattice 
and triangular lattice that report possible realization of 
spin nematic phasesj^"— Among frustrated ferromagnets, 
the spin-i square lattice J1-J2 model with ferromagnetic 
Ji is a prototype minimal model and attracting interest 
recently)^"— ii^"— The model Hamiltonian 
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S-j ■ Sn 
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S-j ■ Sn 



(1) 



consists of nearest-neighbor ferromagnetic exchange 
Ji (< 0) and competing next-nearest-neighbor antifer- 
romagnetic exchange J2 (> 0). When the antiferromag- 
netic coupling is much stronger than that of the ferromag- 
netic coupling (|Ji| ^ 2J2), the ground state exhibits 
a coUinear antiferromagnetic order, (Sj) = {—ly'^m or 
{—lyym with j = [j^^jy)^ While the ground state in 
the opposite limit (|Ji| ^ 2J2) is the fully polarized 
ferromagnetic state. The preceding exact diagonaliza- 
tion studies^rJi suggest the existence of the spin nematic 
phase between these two magnetic ordered phases. 

Motivated by these studies, the present authors re- 
cently formulated a fermionic mean-field theory for quan- 
tum frustrated ferromagnets where they proposed to 
describe the quantum spin nematic phase as a spin-triplet 
variant of spin-rotational symmetric quantum spin liq- 
uids. The state suggested by this mean-field theory is 
a 'mixed' resonating valence bond (RVB) wavefunction, 
where all the half-spins in any lattice points belong to ei- 
ther the singlet valence bonds on antiferromagnetic links 
or the spin-triplet valence bonds introduced on ferromag- 
netic links. The wavefunction is given by a superposi- 
tion of different partitioning of spins into either singlet or 
triplet bonds, such that the wavefunction, on the whole, 
has no preference for any specific valence bond configu- 
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ration. 

Unlike Neel-ordered states or dimer states, this mixed 
RVB state preserves the lattice-translational symmetry 
of the square lattice, indicating its 'quantum spin liq- 
uid' like character. In contrast to the spin-singlet RVB 
state, however, the spin-1 moment in the triplet va- 
lence bond breaks spin-rotational symmetry, quantum- 
mechanically rotating within a specific plane. Thus, un- 
like spin-rotational symmetric quantum spin liquids, this 
mixed RVB state is accompanied by the breakdown of 
global spin-rotational symmetry. In fact, this sponta- 
neous symmetry breaking manifests itself as the ordering 
of the quadrupole moment ^i^ '^^i^^ 

instead of the ordering of the spin dipolc moment, i.e., 
(Sj) = 0. In Eq. j and m denote two adjacent 
lattice sites connected by a ferromagnetic link, and ^ and 
u are the spin indices. Having both quantum spin liquid 
character and also symmetry breaking phase character, 
this mixed RVB state is regarded as the quantum spin 
analogue of 'liquid-crystal' like state of matter 

To investigate the nature of this unconventional quan- 
tum spin state, we study, in this paper, the spin-^ 
square lattice J1-J2 Heisenberg model, using the vari- 
ational Monte Carlo method. Our method is based on 
the fcrmionic representation of spin-i operators Sj^^ = 

y}jc^,.Ufjj itf = l,2,3)^^>25 whcre_4„ denotes 
the fcrmion creation operator at the site j with spin a 
and the Pauli matrices. The representation becomes 
exact, when one and only one fermion is located at every 
site, i.e., /J^/j.a = 1. At the mean-field level, this local 
constraint on fermions' number is replaced by their cou- 
pling with the chemical potential, so that the only the 
global constraint is taken into account. 

In this fcrmionic representation, the exchange interac- 
tion between two localized spins is written as a four-point 
interaction, which the mean-field theory decomposes into 
the pairing fields between two adjacent fermions. The 
antiferromagnetic exchange interaction is decoupled in 
terms of the spin-singlet pairing fields aai^"— 

Si-Sj ~ ^(iXijP + hijP) + const 

+ \{- ^ji -flafj.o. - Vji fla[^'^2]aP.flp + ll-C ) ,(3) 

where Xij denotes the particle-hole (excitonic) pairing 
field, 

Xij = ifljj^c), (4) 
and rjij is the particle-particle (Cooper) pairing field, 

Vij = (/i,a[-icr2]Q/3/j,/j)- (5) 

On the other hand, the competing ferromagnetic ex- 
change interaction should be decoupled into the spin- 



triplet channel in its own right jis 

-Si-Sj ~ ^{\Eij f + \Dijf) + const. 

1 ^ 

+ i E ( - ^J-^'A. fLH'^f^'^^Uflp + h.c.) , (6) 

where Djm = {Djm.i, Djm,2, Djm.s) denotes the 
so-called d-vector of the spin-triplet Cooper pairing 
field^^ 

Dij.f_i = {fi,a[i(^2(^ti\afifj,fi), (7) 

and Ejm = {Ejm,i, Ejm,2, Ejm,^) is the excitonic coun- 
terpart of the d-vcctor, 

Eij,,^{fUc7,Ufj^f,). (8) 

Physically speaking, this decoupling is because the 
ferromagnetic exchange interaction between two spin 
halve o^^'^^ usually prefers the formation of their spin- 
triplet valence bond, \S ~ 1,5'^ = 0} (/i = 1,2,3), in- 
stead of the singlet valence bond. Indeed, the d-vector 
associated with the triplet pairing specifies the direction 
along which the spin triplet state has the zero magneti- 
zation, i.e. 15* = 1, (e^ ■ S) ^ 0) with = Djm/\Djm\- 
One can easily see at the mean-field level that the triplet 
pairing fields in Eqs. ([7]) and ([5]) induce the following 
quadrupole moment)^ 

{Krj)=-l{Eij.,E*ij,,-l5,,\Eijn 

- \{Dij,^,Dlj^„ - i<5^.|Ajf ) + h.c. 

The previous fcrmionic mean-field analysis shows that 
the present square- lattice model, Eq. ([1]), has five dif- 
ferent spin-triplet pairing states as its saddle point solu- 
tions)^ (i) Z2 planar state, (ii) Z2 polar state, (iii) SU (2) 
chiral p-wave state and (iv) 'flat-band' state, where the 
first state is accompanied by a 'coplanar' or 'd-wave' con- 
figuration of the quadrupole moments while the second 
and third ones support 'collinear' configurations of the 
quadrupole moments. Among them, Z2 planar phase and 
SU (2) chiral p-wave phase appear, having the lowest en- 
ergy, in the finite range of competing coupling regime 
between the ferromagnetic phase and tt-Aux phase in the 
mean-field phase diagram. 

In this paper, we investigate the nature and energet- 
ics of the projected BCS wavefunctions constructed from 
these mean-field pairing states. As the local constraint of 
the fermion density is not strictly observed in the mean- 
field theory, the BCS wavefunctions generally range over 
the 'extended' Hilbert space, which allows double occu- 
pancy or vacancy on a single site. To obtain a proper 
trial many-body wavefunction for the spin model, we first 
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FIG. 1: Phase diagram of the spin- 1/2 square lattice Ji~ 
J2 model with ferromagnetic (FM) J\ and antiferromagnetic 
(AF) J2, obtained from variational Monte Carlo simulations. 



project these BCS wavefunctions onto the physical spin 
Hilbert space, imposing the single-ferniion condition on 
every site. By minimizing energies of these projected BCS 
wavefunctions in a variational way, we obtain their opti- 
mal energies. Comparing these energies with those of the 
ferromagnetic state and the collinear antiferromagnetic 
state, we argue that only a projected Z2 planar state 
becomes most energetically favorable among the other 
competing states in a finite range of the intermediate 
coupling regime, 0.42|Ji| < J2 ^ 0.57|Ji|. Our results 
are summarized in Fig. [TJ 

Based on this observation, we further study the char- 
acter of the projected Z2 planar state. Specifically, we 
clarify the irreducible representations of this many-body 
wavefunction under the point group symmetries of the 
square lattice, and argue that this state is actually ac- 
companied by a 'd-wave' ordering of the quadrupole mo- 
ments. All irreducible representations and the d-wave 
character are totally consistent with the nature of the 
spin nematic phase suggested by the previous exact diag- 
onalization study^2. This agreement in combination with 
the energetics suggests that the projected Z2 planar state 
is indeed realized in the intermediate coupling regime of 
the square lattice J1-J2 model. To give a direct physical 
characterization to this quantum spin nematic phase, we 
further calculate the static correlation functions in this 
projected BCS wavefunction in a large system size (100, 
144, and 324 sites). By that, we found that the wavefunc- 
tion exhibits strong antiferromagnetic fluctuation with 
the wave vectors k — (tt, 0) and (0, tt) associated with the 
proximate collinear antiferromagnetic phase, though the 
finite size scaling suggests that, the state does not possess 
any staggered sublattice magnetization in the thermody- 
namic limit. 

The remaining sections of this paper are organized as 
follows: In Sec. II, we briefly review spin-triplet pair- 
ing states obtained by the previous fermionic mean-field 
analysis. We also extend these mean- field solutions into 
the case under finite external Zeeman field. By this ex- 
tension, the projected 'flat-band' state turns out to be 
the fully polarized ferromagnetic state. We also find 
that all the d-vectors in the spin-triplet pairing states 
are lying within a plane perpendicular to the applied 



field. This feature guarantees the 'spin-nematic' char- 
acter of the projected BCS wavefunctions constructed 
from these pairing states. In Sec. Ill, we give a general 
expression for the projected BCS wavefunctions having 
both spin-triplet and spin-singlet pairings. Based on this 
expression, we have optimized numerically the energies 
of (i) projected Z2 planar state, (ii) projected Z2 po- 
lar state, and (iii) projected SU (2) chiral p-wave state. 
In Sec. IV, after briefly explaining the method of opti- 
mization and Monte Carlo simulations, we compare their 
optimized energies with other competing states such as 
the collinear antiferromagnetic state and the fully po- 
larized ferromagnetic state. Sections V and VI contain 
discussions about the nature of the projected BCS wave- 
functions. In Sec. V, we argue that all the projected BCS 
wavefunctions studied in this paper have a spin-nematic 
property, i.e., ordering of the quadrupole moments with- 
out any ordering of spins. We show in particular that 
the projected Z2 planar state is accompanied by a d-wave 
spatial conflguration of ordered quadrupole moments. In 
Sec. VI, we discuss the behavior of the static correlation 
functions of spins and qudrupole moments calculated in 
this projected Z2 planar state. Section VII is devoted to 
the summary and discussion. 

II. MEAN-FIELD ANSATZ UNDER THE FIELD 

The J1-J2 frustrated ferromagnetic square lattice 
model has four types of spin-triplet pairing states as 
the saddle point solutions of the SU{2) fermionic mean- 
field theory: (i) Z2 planar state, (ii) Z2 polar state, (iii) 
SU (2) chiral p-wave state and (iv) 'flat-band' state, all 
of which possess the same translational symmetry as the 
square lattice. We describe in this section how these 
triplet pairing states are deformed under external Zee- 
man field. Wc will sec that all the d-vectors in the states 
(i), (ii) and (iii) are restricted within a plane perpen- 
dicular to the applied magnetic field. Because of this 
arrangement, the mean-field Hamiltonian for these three 
states are invariant under the spin 7r-rotation around the 
field, when combined with the staggered gauge transfor- 
mation /j.Q — > {—iy^~^^yfj^a- This symmetry property 
actually gives the spin-nematic character to the corre- 
sponding projected BCS wavefunctions (see Sec. V). 

A. Z2 planar state 

In the Z2 planar state in the absence of magnetic 
field, the nearest-neighbor ferromagnetic bonds support 
a coplanar configurations of the d- vector, e.g., 

\ DSf,,2 i^ = 3 + ey), ^ ' 

Eij,^ = (9b) 

with Ex = (1,0) and By = (0,1), while the next- 
nearest neighbor antiferromagnetic bonds support the 
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FIG. 2: (Color online) A schematic picture of the spatial 
configuration of d-vectors in (a) Z2 planar state, (b) Z2 polar 
state and (c) SU{2) chiral p-wave state on the square lat- 
tice under the field. The d-vectors in the Cooper channel are 
drawn by (blue) arrows and those in the excitonic channel are 
given by (red) round head arrows, both of which are on the 
nearest-neighbor ferromagnetic bonds. Since all these three 
states preserve the translational symmetry of the square lat- 
tice, we show their configurations in the unit cell. All the 
d-vectors are lying within the plane perpendicular to the ap- 
plied field. These three states are invariant under the spin 
TT-rotation around the field, combined with staggered gauge 
transformation, which guarantees the spin-nematic charac- 
ter of the corresponding projected BCS wavefunctions (see 
Sec. V). 



'staggered-flux' configurations of tlie spin-singlet pair- 
si 



when i = j + Bx i ey. The corresponding Bogoliubov 
de-Gennes (BdG) mean-field Hamiltonian has the form 

^planar = ^ j^"^ (^-^jla ["'s] a/j/J-l-ex ,/3 f j ,afj+ey ,a 
j 

(11) 



J2 

4 



<T = ± 



h.c. 



This Z2 planar state is energetically degenerate under 
spin SU{2) rotation. 

In the presence of magnetic field, spin z-component 
couples to both magnetic field and mean fields from the 
surrounding spins. This effect can be captured by adding 
the term — /lofrX^j to the mean-field Hamiltonian. A 
direct energy optimization suggests that, under the ex- 
ternal Zeeman field, the coplanar plane of the d- vectors 
is restricted within a plane perpendicular to the field. 
Moreover, the magnetic field gives rise to the excitonic 
triplet pairing, whose director vectors are also perpendic- 
ular to the field, replacing Eq. (j9b|) with 



Ei. 



iES 



fi,2 



{i=j + ey). 



(12) 



The mean-field Hamiltonian has the following form under 
the field 



^planar ^ ^planar + ^ ] 



\Jl\ 



The state preserves the translational symmetry of the 
square lattice, so that the mean-field Hamiltonian is 
Fourier-transformed as 



planar 



E fl{ 



\Ji\ 



fc;fe„>0 



■^£^(52:723 + Sy725) " J2XCxCy'r4 
J2'nSxSyf2 + ^/lcff735}/fc, (13) 



{kx,ky), fl — ifl^^i fk^^i f~k,ti f~k,l), fk,a — 

^ .e*''J/j,Q, Sa = sinfca, and Ca = cosfca with a = 



where k 

X, y. The 4x4 7-matrices are defined as 71 = <T2 <Ti, 
72 = 0-200-2, 73 = 0-2(8)0-3, 74 = 0-3(8)0-0, 75 = 0-18)0-0, 
and -fjm = —iljlm, where the 2x2 Pauli matrices in 
front of the 8)-mark is for the particle-hole space, while 
the others are for the spin space, e.g. 



Xij = X, 



(10) 



71 = 





icri 



-ICTi 
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B. Z2 polar state 

The Z2 polar state at zero field takes a collinear con- 
figuration of the d- vectors on ferromagnetic bonds, e.g., 



Eij^fji — 



(14a) 
(14b) 



for i ^ j + Ba (a — x,y), while it supports the same 
staggercd-flux configuration of the singlet pairings as in 

Eq. (nnD. 

In the presence of magnetic field, the external Zeeman 
field restricts the c?-vectors within its transverse direc- 
tions the same as in Eq. (|14ap and brings about the ex- 
citonic spin-triplet pairings, replacing Eq. (jl4bp with 



Ei 



iE5 



>,2 



(15) 



for i = j + Ba {a = x,y). The d- vectors of the induced 
excitonic spin-triplet pairings are perpendicular to both 
the external field and the c?-vectors of the Cooper chan- 
nel. The mean-field Hamiltonian for the Z2 polar state 
has the following form under the field. 



\Jl\ 



-a 

[0'2]a/3/j-|-e„,^i 



a—x.y 



<y=± 



h.C. 



(16) 



The antiferromagnetic exchange interaction supports the 
'uniform-RVB' configuration of the singlet pairingS)^ 



C. SU{2) chiral p-wave state 

In the SU{2) chiral p-wave state, the d- vectors on 
the nearest neighbor x-link and that on the y-link are 
collinear with each other. One of the two acquires an 
additional phase factor i, compared with the other. 



{i=j + By), 



(18a) 
(18b) 



Eqs. P^ -b) hold under any magnetic field, pro- 
vided that the d-vector is perpendicular to the field. 



Xij = X 



(19) 



for i=j + exiey. The mean-field Hamiltonian for the 
SU (2) chiral p-wave state takes the following form. 



E 



p-C'(/j>[^3]a,3/. 



t 



(J=± 



+ h.C 
or equivalently, 

"^chiral p wave 



(20) 



1^1 1 



where ai denotes the uniform temporal gauge field that 
has a finite value in the Z2 polar state (see Ref. [isl ). Or 
equivalently, 

^polar = 5Z fk\- "172 - -^(S^r + Sy){D-f^ - E^25) 
k-ky>Q ^ 

- ■hxCxCy'^i - J2VSxSy'r2 + ^/icff735 l/fc- (17) 



k-.ky>0 

- J2XCxCy74 + ^/ieffTSsj/fc. (21) 



D. Fully polarized state out of 'flat-band' states 

The 'flat-band' states have only spin-triplet pairings 
on the ferromagnetic bonds, while no spin-singlet pairing 
on the antiferromagnetic bonds. According to our pre- 
vious workfi^ this state achieves the best mean-field en- 
ergy among others in the strongly ferromagnetic regime 
(|Ji| » J2). In the absence of the external field, the 
triplet pairings in the 'fiat-band' state are most gener- 
ally characterized by a C/(l) phase, 9, three orthogonal 
unit vectors {ni,n2,'^3} in the spin space, and two or- 
thogonal unit vectors {1711,1712} in the gauge space as 
follows 



cos 9 rii (mi_3 + i m\^2) 
sin^^ni(TO2,3 + «"^2,2) 



(i = i 
(i = 2 



By), 



cos 9 (n2 
sin 9 (n3 



- i riimi^i) 
inim2,i) 



{i=3 
{i=3 



+ ex), 



where 



ni ■ 112 ^ 712 ■ ^ ■ ni = 0, 
mi ■ m2 = 0. 



(22a) 
(22b) 



(23) 
(24) 



The energy dispersion of the Bogoliubov particle com- 
prises two bands, which are totally flat in the momen- 
tum space and energetically separated by 2| Ji| from each 
other. In the remaining part of this section, we will argue 
that this state actually reduces to a fully polarized fer- 
romagnetic state, once an infinitesimally small Zeeman 
field is applied. 

Under the Zeeman field, it is energetically favorable 
that all the d-vectors are perpendicular to the field. To 
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make this compatible with the orthogonahty condition 
Eq. (1231), the C/(l) phase 6 is going to be locked in 6* = |? 
with I G 1. Namely, when 6 = only ris and ni are 
required to be perpendicular to the field, while, in the 
case of = or TT, only rii and n2 are perpendicular to 
the field. Thereby, this locking reduces the 'fiat-band' 
state into a decoupled one-dimensional fermion states 
running along either x-link or y-link. For example, when 
6 — 0, one of the flat-band states under an infinitesi- 
mally small field can be described with ni = (1,0,0), 
n2 = (0,1,0), and mi = (1,0,0). The corresponding 
BdG Hamiltonian is given as 



Hflat = J2 



(25) 



Or equivalently, 



jy = l —Tr<k^<TT 

+ ^/lcff735}/t,.j„, (26) 

^^'^^■'^ fL,jy = (4l j,,T'Alj,,;'/-fexj„.t'/-fcxj„,;) and 



s ^ Et:=i With L^Ly = N. 

When projected into the spin Hilbert space, the 
ground-state wavefunction of Eq. (j26p reduces to a fully 
polarized ferromagnetic state. To see this, note that it 
is given by a composite of decoupled one-dimensional 
fermionic states running along the cc-link; |^fiat) = 



I\jy=i For every jy 

by 

i^.j = n 



1, 



3y/ 



is given 



- cos ^ e ■» + 



sm — e f " 
2 



|o), 



(27) 



where |0) denotes the vacuum state of the fermions. We 
have omitted the index 'jy\ fl^ ^ fl^a^ since the 
argument holds for each jy independently. Note also that 
the quantization axis of the spin was taken along the field. 
The angle cj) is defined as 



tan 



in the range — f < 

The inner-product between \'^jy) and an Ising spin 

configuration \{ajj) = {11^^=1 /I.<t,J|0) is given by a 
determinant of the L r x matrix, 



(KJI*.J=det 



4 e 



4 e 



2+- 



flcr, e 4 e 



,i27r(-l+^) 



i27r(-2+^) 



,i27r(-L^ + - 



where cr = ±1 (t, 4- respectively) with = — (5cr,i cos ^ + 
i5cr,-isin^. This determinant becomes non-zero, if and 
only if all the spins arc pointing upward or pointing 
downward: otherwise, the x matrix always has 
two adjacent column-vectors which are parallel to each 
other in the Lx dimensional space, at a domain wall with 
(cT;, (T;+i) = (—1, 1). Accordingly, we have 

+ (sinf)'^n'5...,-i. (28) 

Note that cos | > sin| > for /leff > 0. This sug- 
gests that, when normalized in the thermodynamic limit. 



Eq. ((27)) always reduces to the ferromagnetic state, 

r^'"^ //,t/ut, J ^..)-lt,t,---,t) (29) 
for any = 1, ■ . • ,Ly. 

III. TRIAL MANY-BODY WAVEFUNCTION 

As shown in the previous section, the projected 'flat- 
band' state under an infinitesimally small field reduces to 
the fully polarized ferromagnetic state. As such, we re- 
gard that the projected 'flat-band' state in the absence of 
the field is the trivial ferromagnetic state, whose energy 
is exactly estimated as — ^(|^i| — J2) (per site). Hence 
we will focus on the character and energetics of the other 
three spin-triplet RVB states, (i) Z2 planar state, (ii) Z2 
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polar state and (iii) 5(7(2) chiral p-wave state. To this 
end, we will construct in this section the projected spin- 
triplet BCS wavefunctions out of their respective mean- 
field Hamiltonians. 



A. Projected spin-triplet BCS wavefunctions 

Let us first derive a BCS 'many-body' wavefunction 
for the BdG Hamiltonian which has both spin-triplet and 
spin-singlct pairings and hopping integrals. Suppose that 
we have a mean-field Hamiltonian TL = fl.Hkfk 
and a 4 X 4 matrix Hk is diagonalized by a unitary trans- 
formation Uk- We typically use Eqs. (|13ll7l2ip for U. 
For these Hamiltonians, the eigenvalues always appear 
in the particle-hole pairwise manner, 



and 



Hk Uk — Uk 



Afe.i 



^fe,2 



-Afc.2 



(30) 



As such, without loss of generality, we can assume Xk,j 
{j = 1,2) to be positive (semi-)definite. Defining the 

Bogoliubov particle 7^^'- {j = 1,2) as 



7L 7I2 7-fc,2 7-fc,i 



= ( /fe,t -fk.l f-k.t f-k,i ) Uk, 



(31) 



we obtain 



^ = X! {^kjli^jlkj - \k,]l-k,jllk^j]- (32) 

j = l k:ky>0 

Since Xkj > 0, the mean-field ground state wavefunc- 
tion |g.s.) is a vacuum of the Bogoliubov particles, i.e., 
7fcj|g.s.) = 7_fc,j|g.s.) = for any k and j, which leads 
to' 



|g.s.) oc Jl {7fc,i7fc,27-fc,27-fc,i}|0). (33) 

k:ky>Q 

Substituting Eq. (PT|) into Eq. , one can easily obtain 
Y[ {7fc,i7fc,27-fc,27-fc,i}|0) 

k;ky>0 

fc;fey>0 

+'^fc/-fc,f^fc-T + ^'kf-k,ifk.-\ + ^kf^_k.\f-k,ifk.-\fk,i]\'^) 

(34) 



with 



n { - [Uk\iAUk\2A + [UkhAUk^A) 



ky>0 



{[Uk]lAUk]l2-[Uk]l2[Uk]l,} 



(35) 



[Uk]UUk]h-[Uk]l2[Uk]h 

[Uk]uuk]i2-mi2[Uk]h' 

[Uk]*l,l[Uk]l,2 - [^fc]l,2[^fc]4,l 
[Uk]*i,i[Uk]2,2 - [Uk]*i^2[Uk]2,l' 
[Uk]2A[Uk]*3,2 ~ [Uk]2,2[Uk]t,l 



Ck = 



[Uk 


li[Uk 


2,2 ^ 


[Uk]l2[Uk 


2,1 


[Uk 


2,1 [^'^ 


4,2 ^ 


[Uk]i2[Uk 


1.1 


[Uk 


li[Uk 


2,2 ^ 


[Uk]l2[Uk 


2,1 


[Uk 




4,2 ^ 


[Uk]l2[Uk 


4,1 


[Uk 


li[Uk 


2,2 ^ 


[Uk]l2[Uk 


2,1 



(36) 



Notice that Ck = a/c^'fc — '^'k^k- This makes it possible to 
exponentiate the right hand side of Eq. ([34]) as 



|g.s.)EEexp[ ^ [tk]^^flk,jl,] |0) 



with 



tk — 



b'k bk 



(37) 



(38) 



Equation (|37p generally has a finite weight not only 
on physical (i.e. spin) Hilbert space but also on those 
fermionic states having either double occupancy on a sin- 
gle site or an empty site. To obtain a variational many- 
body wavefunction in the physical spin Hilbert space, we 
need to project out these unphysical fermionic states, im- 
posing 'single-fcrmion condition' on every site; 



(39) 



where V stands for the projection operator onto the phys- 
ical spin Hilbert space. The projected BCS wavefunction 
l^'ct) depends on the pairing and hopping fields encoded 
in the BdG Hamiltonian, such as D, E, V ^-^id /loff. 
The characteristic of the mean fields is symbolically rep- 
resented by the subscript a. 

Within the spin Hilbert space, the wavefunction is ex- 
pressed by its inner-product with an Ising spin configu- 
ration, KcTj}) = {Ilj /j o-j This product generally 
reduces to a Pfaffian,i22i^i2^ 



({a,-}|vI/„)=Pf[X„({a,-})], 



(40) 



where Xa{{(Jj}) denotes the N x N antisymmetric ma- 
trix given by^i^ 



k;ky>0 



Note that the boundary condition for the momentum k 
remains arbitrary in Eq. (|4ip . To fix this arbitrariness, let 
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us require that the spin wavefunction given by Eq. pO 
is an eigenstate of translations, 



(W„(,-)}|*a)=e''=(K-}|*.), 



(42) 



where Ta denotes the lattice translational operation by 
Ba, i.e., Ta{j) = j + Ba (a = x,y). In fact, the pre- 
ceding exact diagonalization studies^ suggest that the 
states with non-zero Q vectors are unlikely realized in 
any intermediate coupling regime of the present Ji- 
J2 model, so that we impose the translational invari- 
ance (e*^-,e'^f) = (1,1) on Eq. (gH). To satisfy this 
requirement, the fermion's momenta in Eq. (j4ip have 
only to observe either the anti-periodic boundary con- 
dition (APBC), i.e., ka = {2na — l)7r/La with Ua = 
—La/2 + 1,.-. ,La/2, or the periodic boundary condi- 
tion (PBC), i.e., ka = 2naTT / La- For the two-dimensional 
models, the trial wavefunctions have four options. When 
both kx and ky satisfy the anti-periodic boundary condi- 
tion, the total momentum carried by our trial spin wave- 
function is indeed at the F-point, i.e. (e'^'", e'^") = (1, 1); 
to see this, one has only to relate Xa{{crT^(^j)}) with 
X.a{{aj}) in terms of a certain elementary row/column 
operation Oa, 

where Oa exchanges site indices of Xa{{<7j}) according 
to the lattice translation operator Ta- Similarly, when 
kx satisfies the periodic boundary condition while ky 
does the anti-periodic boundary condition or vice versa, 
the momenta carried by the projected BCS wavefunc- 
tion can be shown to be (e'^=", e'^") = ((— 1)^^""^)^^ , 1) 
or (e*^"' , e*^" ) = (1, (— 1)(^==~^)^!'), respectively. In what 
follows, we only consider the systems with even length 
Lx and Lj,, to impose the translational invariance of the 
total wavefunction, i.e., (e*^^,e*^f) = (1, 1). 

Notice also that, when both kx and ky observe the 
periodic boundary condition, any of the projected BCS 
wavefunctions derived from Eqs. (|13ll7l2ip cannot be ex- 
pressed in terms of a single Pfaffian. This is roughly 
because, being either d-wave or p-wave, all the pairing 
fields in the Cooper channel always vanish at the four 
time-reversal invariant momentums points (0,0), (0, tt), 
(tt, 0) and (7r,7r), where a 4 x 4 BdG Hamiltonian re- 
duces to a 2 X 2 Bloch Hamiltonian having no anomalous 
part. As a result, the state-basis representation of the 
projected BCS wavefunction becomes relatively cumber- 
some. In this paper, we study only those projected BCS 
wavefunctions derived based on the other three bound- 
ary conditions. Following the standard literaturejii^ we 
name the projected BCS wavefunction defined with the 
APBC in the both direction as the 'wavefunction in the 
(tt, 7r)-topological sector' and that with the PBC in one 
direction and the APBC in the other as the 'wavefunction 
in the (0, tt) or (tt, 0)-topological sector.' 



B. Quantum spin number projection 

Our Hamiltonian has the global SU (2) spin rotational 
symmetry, while the trial wavefunctions constructed from 
spin-triplet pairing states explicitly break this continuous 
symmetry by hand. Such a symmetry breaking is sup- 
posed to occur only in the thermodynamic limit. The 
ground-state wavefunction in a finite-size system can be 
always identified as an eigenstate of the symmetry groups 
of the Hamiltonian. Accordingly, it is naturally expected 
that the energy of the trial state will be further improved, 
when the state being projected onto the eigenspace of an 
appropriate quantum number associated with the spin- 
rotational symmetries. Thus, we also consider as our 
trial state the projections of the spin-triplet BCS wave- 
functions with the quantum spin numbers. 

The spin projection operator which filters out a state 
with the total spin S = L and the z-component of the 
total spin 5'^ = M has a forn^-^ 



Vs^=m'Ps=l 



2L + 1 



2ir 



da / dp sin /3 / dj 



2tt 



X Fl(cos/3) e'"(^^-*^)e'^^''e*^^% (43) 



where Pl denotes the L-th Legendre polynomial, and 
T's=L and 'Ps^=m denote the projection operators filter- 
ing out a state with the total spin S = L and a state with 
the z-component of the total spin Sz = M , respectively. 

Combining this with Eqs. ([37 |) -(|4l ]) . we obtain the pro- 
jected BCS wavefunction with the quantum spin number 
projection as 



{{<Jj}\rS^=MVs=L\^c 



2L+1 
47r 



d/3sin/3 / d7Pi(cos/3) Pf[X„({a,-};/3,7)] (44) 



under the condition i CTj = M. Here, the N x N 
antisymmetric matrix Xa{{<yj}; /3, 7) is defined the same 
as in Eq. (|41l) with the 2x2 matrix t{j,l) being redefined 
in a rotated spin frame; 

[X^i{aj};Pn)],,HVp^,tijJ)V,:^]^^^., 

-[y,,,*(i,j)F^;]^^ (45) 



where 



cos ^ — sin ^ 
sin ^ cos ^ 



e'i 
e-^i 



To integrate over (3 and 7 numerically in Eq. (|44|) . we 
employ the Gauss-Legendre quadrature. When project- 
ing into the singlet space, i.e. 5' = 0, with the system 
size = 6 X 6 ~ 12 X 12, we typically used 10 ^ 16 
mesh points for the integration over j3 and 10 ^ 20 mesh 
points for that of 74S 
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IV. ENERGY OPTIMIZATION AND 
ENERGETICS 

In this section, we optimize the energies of the pro- 
jected (i) Zi planar state, (ii) Zi polar state and (iii) 
SU(2) chiral p-wave state and compare their minimized 
energies with those of the ferromagnetic state and the 
coUinear antiferromagnetic state. Specifically, we have 
numerically calculated the expectation values of the en- 
ergy for these projected BCS wavefunctions, taking the 
quantum spin number projection onto the subspace with 
either S*^ = or S" = 0; 



E. 



(*«|7's.=o|*«> 



and 



E: 



(*«|7's=o|«'a) 



(46) 



(47) 



We have further optimized these energies, tuning the 
variational parameters a encoded in the original BCS 
wavefunctions, such as D, E, x, 77 and h^ff. For this 
optimization, we have employed the so-called stochastic 
reconfiguration methodi^i^ 



small positive value A. Meanwhile, the SR method de- 
termines the optimal direction, by minimizing the en- 
ergy on the contour- (super) sphere of the equal quantum 
distance. A variational principle with constraint dic- 
tates that the optimal direction thus defined is given by 
5aj = ^J2m\-'^a^h,m'^oi^-^of It IS empirically recognized 
that the modification in terms of the metric tensor [Sa] 
substantially improves the optimization efficiency, espe- 
cially when the tensor has a highly non-flat structure in 
the variational parameter space<^2^ 

The numerical evaluation of the metric tensor and the 
gradient vector requires the summation over all the Ising 
spin configurations in the physical Hilbert space, 

[Sa]ni,n ^^Yl (^m,{<T,- } C",{<t,- } +C.C.)w{<,.} 

- ^ReOjn.{a^}W{a.} ^ ReO„,{^^} , (51) 



2 RcO, 



where 



A. Stochastic reconfiguration method 

Here we briefly review the stochastic reconfiguration 
(SR) methodi^i^i In this optimization method, a usual 
steepest descent (SD) method is modified in such a 
way that information of the 'quantum distance' between 
wavefunctions is included. The quantum distance is cho- 
sen to be the square distance between two normalized 
wavefunctions defined in two different parameter points, 
say a and a -\- Sex, in the form 



^SR — 



with 



I - l^-a) (48) 



Regarding \Sa\ as a small quantities, we can expand this 
quantum distance in terms of Sa, 



^SR 



ajSa,, 



(49) 



where the metric tensor [Sa] is defined in the variational 
parameter space as 



C.C.. 



(50) 



In the standard steepest descent (SD) method, the vari- 
ational parameters are changed along the gradient of an 
energy, 5a j = XdajEa with E^ = (5'a|-f^|^a) and a 



O. 



l({^j}|^c.)P 





H\Wj]) 







(53) 
(54) 

(55) 



The SR method replaces this extensive summation by 
the statistical average where w^^^-^ is regarded as a prob- 
ability density of the corresponding statistical ensemble. 
Specifically, we numerically create a Markov chain in 
which a binary configuration {ctj} is statistically gen- 
erated with the probability w^^^-^. In the statistical en- 
semble thus defined, observables defined in Eqs. (|54j|55p 
are numerically evaluated; 



1 



[S, 



01 



c.c.) -i(0„+c.c.)(0„ + c.c.). 



c.c.)(C'™ +C.C. 



To obtain the metric tensors and gradient vector, we usu- 
ally take 1000 ^ 4000 samplings per site. For a set of 
optimized variational parameters, we evaluate the en- 
ergy (sec. IVB) and the correlation function (sec. VI), 
where we typically use 10^ samplings. As for the pro- 
jected BCS wavefunction with the quantum spin number 
projection, one has only to replace \^cx) in Eqs. (|53|) - (|55)) 
by Vs^=m'Ps=l\^ol)- Those who are interested in the ac- 
tual evaluation of these observables can consult Refs. [sgI 
and Hi. 
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FIG. 3: (Color online) Energy comparison in the square 
lattice J1-J2 model as a function of Q with (Ji,J2) = 
(— sin cos 0). (a): Optimized energies of the planar states 

[(A) Ps = o|*planar> aud (B) Ps, =0 | *planar) ] , the polar StatC 

[(C) Ps,=o|vE'poiar)], the TT-flux State [(D) Ps,=o|*x-tiux)], the 
p-wave chiral state [(E) Psj=o|*I'chirai)], and the collinear an- 
tiferromagnetic state (F). The exact energies of the ferro- 
magnetic state (G) and the isolated dimer state (H) are also 
shown. The projected state Ps=o|*I'pianar) is calculated in 8x8 
spin system and Ps^=o|^' planar) is in 10 x 10 spin system, (b): 
A part of figure (a) is enlarged, (inset): Optimized param- 
eter values of D, X) and r] in 7^3=0 jvl'pianar) as a function of 
d. These values are rescaled, such that \J D'^ + + = 1- 
The optimal values of the excitonic spin-triplet pairing field 
E and the effective Zeeman field h^s are negligibly small. 



B. Energetics 

Figure [3] shows the optimized energies of the 
projected Z2 planar states (both "Ps-^o I ^planar) 
and 7^s=o I ^planar)), the projected Z2 polar state 
7^s.=o|^poiar)! and the projected SU{2) chiral p-wave 
state 'Ps^=o|^chirai). We also compare these optimized 
energies with the exact energy of the fully polarized 
ferromagnetic state Sforro = — 0.5(|Ji| — J2), and 
the variationally optimized energies of the collinear 
antiferromagnctic stato^ Ecaf — — O.6682J2 and the 
decoupled double tt-Aux stated ^-x tiux = — 0.64J2. The 
collinear antiferromagnctic state we considered in this 



0.11 



0.1 



2 0-09 

I 0.08 



0.07 - 
0.06 



B 



1.1 

[rad] 



1.2 



FIG. 4: (Color online) Comparison between the optimized 
energies of the projected planar states (A) Ps=o|*I'pianar) and 
(B) Psj^ol^pianar), and the ground state energy ^Eed ob- 
tained by exact diagonalization, in the square lattice J1-J2 
model, as a function oi 9 = tan^^ (— Ji/ J2). The system size 
is 6 X 6. 



paper is composed by two 'decoupled' Necl-ordered 
states 



I* 



CAF 



) |1'Ncci)a|1'Ncci) 



IB, 



(56) 



each of which is defined on a non-frustrated square 
sublattice coupled with J2 bonds, say A-sublattice or 
B-sublattice. For I^'nogi), we employed the varia- 
tional wavefunction numerically derived by Liang et.al.;^ 
which was maximally optimized on the antiferromagnctic 
square lattice in terms of the staggered magnetic mo- 
ment, singlet pairing fields, and the Jastrow factor. Note 
that the expectation value of the ferromagnetic exchange 
interaction always vanishes in Eq. (j56|) . though the spin- 
triplet pairing fields connecting the two sublattices could 
possibly decrease the energy in general. In spite of this, 
however, the energy of Eq. (|56|) achieves about 94.5% 
of the exact ground state energy in 36 spin cluster at 
|>^i| = J2, whereas it achieves about 96.3% at Ji = 0. 
We thus regard that, even in the presence of consider- 
able Ji , Eq. ([55]) still gives an appropriate energetics for 
the collinear antiferromagnctic state of the J1-J2 model. 

The energy comparison in Fig. |3] shows that the pro- 
jected Z2 planar state has the lowest energy in a finite 
range of the intermediate coupling regime, 0.417|Ji| < 
J2 ^ 0.57|Ji|, whereas the ferromagnetic state is the 
most stable in the strong Ji regime, J2 < 0.417|Ji|, and 
the collinear antiferromagnctic state is in the strong J2 
regime, 0.57| Jij < J2- The optimal energy of the planar 
state projected onto the Sz = sector in the 6x6 system 
achieves roughly 92% ^ 89% of the exact ground state 
energy obtained by the numerical diagonalization with 
the same system size. The energy becomes further de- 
creased by 2% ^ 3%, when the wavefunction is projected 
onto the S' = space (see Fig. H]) . 

Figure [3] suggests that, contrary to the mean-field 
analysis, the projected chiral p-wave state hardly real- 
izes in any of the intermediate coupling regime of the 
J1-J2 model, at least when the system size is 2n x 2n 
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FIG. 5: (Color online) Size dependence of the optimized ener- 
gies per site for the projected planar states in the LxL square 
lattice J1-J2 model with J2 — 0.55| Ji | (a: black dashed lines), 
J2 = 0.5|Ji| (a: red dashed lines), and J2 = 0.45|Ji| (b). The 
energy unit is taken to be y>^f+J|. The projected states 
are taken as ■Ps^=o|^'pianar) (L = 6 ~ 14) for the data (B, D, 
F, H) and as Ps=o|*pianar> {L = 6,8) for the data (A, C, E, 
G). The horizontal axis is taken as L~^. For Psj=o|*I'pianar) 
with L — 14, we used the same variational parameters as 

those of Ps^^ol^planar) with L = 12. For "Ps^ol^planar) with 

L = 8, we used the same variational parameters as those of 
7'5_o[*I'pianar) wlth L — 6. The data G and H are calcu- 
lated in the (0, 7r)-topological sector and the others are in the 
(tt, 7r)-topological sector. 



(71 = 3,4, 5---). Moreover, the estniiated energy of 
■Ps^=o|^'chirai) IS already 30% higher than those of the 
ferromagnetic state and the projected Z2 planar state, 
so that the situation is unlikely reversed, even when the 
wavefunction is further projected into the singlet space, 
i.e. 'Ps=o|^chirai)- Figurc [3] also indicates that the Z2 po- 
lar state is almost energetically degenerate with the dou- 
ble TT-flux state, which indicates that the spin-triplet pair- 
ing does not lower the energy efficiently in this state. In- 
deed, we observed that the optimized value of the triplet 
pairing field in the polar state is less than 10% of the root 
square sum of all the variational parameters. 

We also show the system size dependence of the opti- 
mized energy per site for the projected planar states in 
Fig. [5] The data for the state 7^5^=0 1 ^planar) indicates 
that the energy has a finite-size correction in the form 
E{L)/L'^ = eg — c/L^. This observation is consistent 
with the existence of gapless Goldstone modes in the Z2 
planar state. Using the random phase approximation, we 
can show that the low energy dispersions of these gapless 
modes are always linear in the momentum;^ which sug- 
gests that the finite-size correction to the ground state 



energy per site decays in the form —c/L^, the same as 
that of the two-dimensional antiferromagnetic Heisen- 
berg model4i 



V. d-WAVE SPIN-NEMATIC CHARACTER OF 
PROJECTED Z2 PLANAR STATES 

In this section, we show that all the projected spin- 
triplet RVB states derived from Eqs. (jl3ll7l2ip gen- 
erally have the 'spin-nematic' properties; ordering of 
quadrupolc moments without spontaneous ordering of 
magnetic dipole moments. In particular, we argue that 
the projected Z2 planar state has a 'd-wave' spin-nematic 
character, or an 'd-wave' quadrupolar order, which is con- 
sistent with the nature of the spin nematic phase sug- 
gested by the exact diagonalization study^. All symme- 
tries of the Anderson tower of spin nematic states are 
clarified from the decomposition of the Z2 planar state. 

To see the 'spin-nematic' character, notice first that 
all the mean-field states discussed in Sec. [Tl] are invari- 
ant under the spin 7r-rotation about the 3-axis. Namely, 
spin-triplet d-vectors in these states are always lying in 
a plane perpendicular to the field (see Fig. 1), so that 
the spin 7r-rotation around the field changes the sign 
of the triplet pairing fields on the nearest-neighbor fer- 
romagnetic bonds, while leaves intact the singlet pair- 
ing fields on the next-nearest-neighbor antiferromagnetic 
bonds. This sign change can be readily set off by the 
staggered gauge transformation /j — >■ /j (— !)•'==+•'!' . The 
whole unitary transformation is expressed as 



U =exp 



X exp 



(57) 



Since the mean-field Hamiltonian is invariant under this 
transformation, the transformed state J7|g.s.) is energeti- 
cally degenerate with the original ground state |g.s.). On 
the one hand, being a vacuum state of the Bogoliubov 
particle, the ground state should be unique, which leads 
to U\g.s.) = e'^lg.s.). Thereby, the projected BCS wave- 
function generally satisfies the following relation 

= {{aj}\V\g.s.) = e-^'{{aj}\VU\g.s.). (58) 

Since the unitary operator U commutes with the pro- 
jection V, the right hand side can be further written as 
follows with the eigenvalues Sz = ^ j "'j i 

({a,-}|*„) = e-''e"«=({a,-}|vI/„). (59) 

One can easily fix the U{1) phase factor, evaluating the 
product between the projected BCS wavefunctions and a 
fully polarized state \{(Tj = 1}) s {Y[j /j^}|0), 

\{Wj = mc.)?= n l^'^l'' 

fc;fej,>0 
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where a'^ is defined in Eq. (|36l) . Supposing that aj, is non- 
vanishing at any discretized momentum point at ky > 0, 
which actually holds true for the Z2 planar state in the 
zero field case, i.e. Eq. PT|) . the projected BCS wave- 
function has a finite weight in the eigenspace of 5*^ = , 
\{{<Tj = 1}|\I>„)| ^ 0. To make this observation compati- 
ble with Eq. ([5^ . the phase factor must take a form, 



(60) 



In general, we can prove Eqs. ()59ll60p more directly, only 
by imposing the spin-7r rotational symmetry onto the 
eigenvectors of a given Bogoliubov Hamiltonian. 

Equation (|59p guarantees the 'spin-ncmatic' character 
of the spin-triplet RVB states. Namely, when combined 
with Eq. (j60p . this identity requires that the wavefunc- 
tions have a finite weight only in the subspace with an 
even-integer Sz for iV = 4Z (Z = 1, 2, • • • ) spin systems, 
whereas only in the subspace with an odd-integer Sz for 
N = 41 + 2 {I = 0,1, ■ ■ ■) spin systems. Thus, the trans- 
verse local magnetization always vanishes in these pro- 
jected spin-triplet RVB states. 



(61) 



while the spin quadrupolc moments in the transverse 
plane are allowed to have a finite value. 



(62) 



where Sj +Sm. + relates to the spin nematic operators in 
the formi ^,-,+5^,+ = (i^|^ - K]^J + 2zK]^^. 

As for the Z2 planar state derived from Eqs. (|llll3p . 
this quadrupole moments obey the d-wave spatial config- 
uration. 



f{R^{j-m))^-f{j-m), 



(63) 



where Rg denotes the space 0-rotation around the axis 
perpendicular to the square-lattice plane. This d-wave 
nature comes from the fact that the planar state is invari- 
ant under the space -l-rotation accompanied by the spin 
^-rotation around the 3-axis (around the field) and the 

gauge transformation /J o- ~^ *(~1)"''"^"'''/J o-- Namely, 
the state is invariant under the following unitary trans- 
formation, 



U' =exp 



X exp 



(64) 



Utilizing this symmetry in the same way as we did for 
the spin 7r-rotation above, one can derive the following 
identity for the projected Z2 planar state, 



TABLE I: Indices of the projected Z2 planar states under 
the point group symmetries of the square lattice, where A'^ 
(multiples of 4) denotes the total number of the lattice points. 



point group symmetries 



"P5^=2n|^planar) 



time-reversal (for n = and zero-field case) 1 

■j-rotation within the lattice (~1)^ 

mirror with respect to the a;-link 1 

translations 1 



To obtain the d-wave character from this identity, expand 
the projected planar state into each S'^-subspace. Under 
Eq. ([5^ . it takes a form, 

l^planar) = ' ' ' + ^^5= =-2 | ^pl anar anar/ 

(66) 

for N = 41 {I = 1,2, ■ ■ ■) spins. Among this set of states, 
the quadrupole operator connects only those two states 
whose Sz differ by 2. 

f{j-m) = J2 

n 

(^planar I 'Ps^=2ri-|-2 Sj^ + Sm,+ "Ps^ =2ri | 'fplanar) • 

From Eq. (pS]) . one of these two states is always even 
under the space ^-rotation, whereas the other is odd. 
This clearly assigns the d-wave spatial configuration of 
the quadrupole moments, i.e. Eq. ([63)) . 

The arguments so far also suggest how to construct a 
trial wavefunction for the so-called Anderson's tower of 
states (or quasi-degenerate joint states) of the symmetry 
breaking d-wave spin nematic order under the field. For 
N = 41 spin clusters, 7^s,=2n|5'pianar) mimic the quasi- 
degenerate joint states (QDJS) of the spin nematic or- 
dered phase, whose representation under the point group 
symmetry operators are listed in Table. I. These repre- 
sentations are actually consistent with those of the Bose- 
Einstein condensate phase of a two-magnon bound state 
near the saturation fieldi^ 

The projected Z2 planar state in the zero-field case is 
time-reversal symmetric; the wavefunction derived from 
Eq. (jlip preserves the following symmetry property, 

Pf [Xpianar({a'n})]* 

= (-1)^ { \{i-^r' } Pf [^Planar({-CT„})] ■ (67) 
3 

To see this relation, notice first that the time-reversal 
operation changes the sign of the triplet Cooper pairing 
fields in Eq. (|TT|) . while does not affect the singlet pairing 
fields. Such a sign change can be readily compensated 



by the previous staggered gauge transformation, f - — >■ 

/j (-1)J-+J«, which imposes the following relation onto 
the BCS gap functions: 



({^j}|*planar) = (- 1 ) " ({(Tfl (j) } | ^-planar) • (65) 



[tfc] — (T2 [<_fe+(7r,7r)]c2- 
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FIG. 6: (Color online) Spin correlation functions of the 
transverse component C±{j) (A and C) and the longitu- 
dinal one C'zz{j) (B and D) in the projected Z2 planar 
state 'Psj=o|*I'pianar) aloug (a) the x-direction (1,0) and (b) 
the diagonal direction (1,1). [A: C±(j,0), B : C^^(j,0), C: 
C±{j,j) and D: Czz{j,j)-] The projected planar state was 
obtained in the J1-J2 model with Ji — —1 and J2 = 0.45 in 
18 X 18 spin cluster, which takes the variational parameters 
as {D,x,v) = (0.40,0.57,0.72). The error-bar is smaller than 
the symbol. 



Or equivalently, 

[-X'planar({o'Ti})] 



x(-l)-^ (-l)-- [Xpia„a,({-<T„})] 



jm 



Noting that Pf [O^AO] = detOPf[A], one immediately 
obtain Eq. (|67p . This equation especially means that, 
as far as l^pianai) is constructed from Eq. ((TT|) . both 
■Ps,=o I* planar) and 7^5=0 1 *pianar) are evcH undcr the 
time reversal operation for N = 41 spin systems, which is 
also consistent with the nature of the spin nematic phase 
suggested by the previous exact diagonalization studiesi^ 



VI. STATIC CORRELATION FUNCTIONS 

Based on energy comparison and symmetry arguments, 
we have argued so far that the projected Z2 planar state 
is likely to be realized in the square lattice 5 = 1/2 
J1-J2 frustrated ferromagnetic model in the interme- 
diate coupling range 0.417|Ji| < J2 < 0.57|Ji|. To 
give a direct physical characterization to this interme- 
diate phase, we discuss in this section the static correla- 
tion functions calculated with respect to Ps^=o I ^planar) 
and Ps=o I ^' planar)- From the energetics, it is clear that 
7's=o|^'pianar) is closer to the symmetric ground state of 



finite spin clusters than the other in the intermediate 
phase. On the other hand, the correlation function of 
7's^=o|^pianar) offcrs a feature of spin-rotational symme- 
try broken spin nematic state. 

Let us begin with the spin correlation function calcu- 
lated with respect to 7^s^=o|^'pianar)- In Fig. [51 we show 
the characteristic behavior of the transverse component 
of the spin correlation function (labeled as 'A' and 'C'), 

C+_(j -m) = i(*planar|^5.=0 {Sj. + Sm.- 

+ Sj-Sm,+ } 'Ps^=o|^'planar), 

and the longitudinal one (labeled as 'B' and 'D'), 

Czzij - m) = ('I'planar|7's^=0 Sj,zSm,z ^S,=0 | * planar ) ■ 

Observing them, notice first that the transverse compo- 
nent of the spin in the A-sublattice in which jx+jy is even 
has no correlation at all with those in the _B-sublattice in 
which jx + jy is odd. More generally, this feature holds 
true for any (projected) Z2 planar state derived from 
Eq. dni) or dni), i.e., 

(^planar| {S'j, + 5'm,- + S'j,_S'm,+ } | ^planar) = (68) 



for y j E A and \/ m E B. Equation (|68p can be un- 
derstood from the symmetry argument. Suppose that 
the spin 0-rotation around the z-axis is applied onto 
all the spins in the A-sublattice, while the spin —6- 
rotation is in the i?-sublattice. The mean-field Hamil- 
tonian for the Z2 planar state is invariant under this 
continuous transformation, so that the staggered mag- 
netization Sa,z - Sb,z = \ Sje-A '^d ~ \ ^j&B '^j is a 
conserved quantity. Applying this staggerecl spin rota- 
tion to the projected Z2 planar states, we obtain 



9iSA..-SB..) 



({''■j} I ^planar) ^ 

for any 6. Equivalently, we have 

({f^j}|*planar) = Ssa^,,Sb.. 9i{<^j}) 



({f^jll^planar) (69) 



(70) 



which immediately leads to Eq. 

Equation ([68)) suggests that, as for the transverse com- 
ponent of the spin correlation function, the next-nearest- 
neighbor antiferromagnetic interaction dominates over 
the competing nearest-neighbor ferromagnetic interac- 
tion. In fact, the transverse spin exhibits a strong antifer- 
romagnetic correlation on each of the 'unfrustrated' sub- 
lattice (square lattice with J2 bonds) , where the function 
decays nearly in a power-law. (see 'A' and 'C in Fig. [5]). 
On the one hand, the longitudinal component always has 
a ferromagnetic correlation between the nearest neighbor 
spins and an antiferromagnetic correlation between the 
2nd neighbor spins. Thus, the spin-frustration among 
the ferromagnetic bonds and antiferromagnetic ones ef- 
fectively suppresses the overall amplitude of the longi- 
tudinal correlation function (see 'B' and 'D' in Fig. [6]). 
Indeed; Czzij) decays quite rapidly and falls below 10~^ 
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FIG. 7: Spin correlation function C{j) in the spin-singlet 
state Ps^ol^'pianar) along the a;-direction (a) and along the 
diagonal direction (b); [(a) C{j,0) and (b) C{j,j)]. The pro- 
jected planar state was obtained for the J1-J2 model with 
J2 — 0.45|Ji| in 12 x 12 spin cluster, which takes the varia- 
tional parameters as {D,x,v) = (0.45,0.55,0.70). 
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FIG. 8: (Color online) A log-log plot of the correlation func- 
tion C{j) in the spin-singlet state 7^5=0 l^'pianar) with the 
same parameter set as used in Fig. [T] We employ only those 
points with ma,x{jx,jy) < 6 and jx +jy = even. The red solid 
line indicates the slope of C{j) ~ ( — 1)"''" |j |~°'^- 



when the spins are spatially separated by more than two 
sites, IjI > 3. 

When the wavefunction is projected onto the spin sin- 
glet space, the static spin correlation function 

Cij -m) = (M/p, 

anar \Vs =0 Sj ■ Sm 

becomes spin-rotational invariant, which 'interpolate' be- 
tween C-i {j) and Czz{j) described above. Namely, as 



FIG. 9: (Color online) (a) Static spin structure factor S{k) 
calculated from 'Ps=o|^'pianar) with the same parameter set 
as used in Fig. [T] The momentum vector fc ranges over the 
1st Brillouin zone, [— 7r,7r] x [— 7r,7r]. (b) Finite size scaling of 
[Afi(7r, O)]'^ = jv+2 ^(''^y 0) ^ function of the linear dimen- 
sion of the system size, indicating that A/L(7r,0) converges to 
zero in the thermodynamic limit within the statistical error 
of the Monte Carlo estimation. 

shown in Fig. [3 the correlation between the spins in 
the A-sublattice and those in the i3-sublattice are ei- 
ther extremely short-range or almost quenched, while 
the correlation within each sublattice exhibits an antifer- 
romagnetic 'quasi-long-ranged' power law decay, which 
is fitted as C{j) - ("IpHjl"'' with 77 = 0.9 ~ 1.0 
for jx + jy = even (see Fig. |5]). Correspondingly, the 
static spin structure factor calculated with respect to 

'7^5=0 I* planar) 

^(fc) = ^e*^-C(j) (72) 
j 

has characteristic peaks at (tt, 0) and (0,7r), while it lose 
its weight at (0, 0) and (tt, tt) (see Fig.[9ja)). This behav- 
ior resembles the structure factor in the coUinear antifer- 
romagnetic phase and seems to be consistent with a re- 
cent exact diagonalization study up to 40 sitesi^ In spite 
of the prominent antiferromagnetic fluctuation at (0,7r) 
and (7r,0), however, the standard finite size scaling fit- 
ting of the static spin structure factor— 1^ suggests that 
the projected Z2 planar state does not have any finite 
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factor of the quadrupole moments 



{n,n) 



(0,0) 



{n,n) 




FIG. 10: (Color online) (a) Static structure factor Tab{k) of 
the quadrupole moments calculated in the projected Z2 pla- 
nar state ■Ps=o|*I'pianar), with the same parameter set as used 
in Fig. [3 for the J1-J2 model with J2 = 0.45| Ji| in 10 x 10 
spin cluster. The (red) solid line depicts the diagonal com- 
ponent Txx{k), the (blue) dotted line the off-diagonal com- 
ponent with minus sign ~Txy{k), and the (black) dashed line 
the other diagonal component Tyy{k). Inset: The momen- 
tum fe is taken along the high-symmetric momentum points, 
(b) Finite size scalings of the long-range order of quadrupole 
moments -^Tab{0,0) as a function of 1/L. The off-diagonal 
component T^y (0, 0) is shown multiplied by minus sign. 



sublattice magnetization M{tt, 0) in the thermodynamic 
limit [see Fig.EIb)]; 

[M{7T, 0)]2 = hm ^^3) e'-'^ ^ 0. (73) 

3 

The correlation functions of the quadrupole moments 
are calculated in terms of the projected Z2 planar state 
in the spin-singlet subspace, 



X A^Zm+e.^5=0l*planar) (74) 

with a,b = X, y. We found that the diagonal components 
thus calculated, D^.j.(j) and Dyy{j), are positive-definite 
for any j, while the off-diagonal component always takes 



r,,(fc)=^e*^'i?afc(j), 



(75) 



both the diagonal components and (minus of) the off- 
diagonal component exhibit prominent peaks at the T- 
point [see Fig. iroTa)]. indicating the d-wave ordering 
character of the quadrupole moments in the projected 
Z2 planar state. In fact, finite size scalings of their peak 
values suggest that the state is indeed accompanied by 
finite quadrupole moments in the thermodynamic limit 
[see Fig.[IOi;b)]. 



VII. SUMMARY AND DISCUSSION 

In this paper, we have investigated the phase diagram 
and the nature of a quantum spin nematic phase in the 
spin-i quantum frustrated J1-J2 model with ferromag- 
netic Ji on the square lattice, describing the ground state 
wavefunction in terms of a spin-triplet pairing state of 
the spinon fields. Our theory is based on the previ- 
ous fermionic mean-field analysis?^ which proposed four 
types of the mean-field solutions. These solutions include 
(i) Z2 planar state, (ii) Z2 polar state, (iii) SU{2) chiral 
p-wave state and (iv) 'flat-band' state, all of which are 
characterized by different kinds of spin-triplet pairings 
of the spinon fields introduced on ferromagnetic bonds. 
Like in usual 'projective description' of symmetric quan- 
tum spin liquids fii^'^i^'^ we construct projected BCS 
wavefunctions out of these triplet pairing states. Per- 
forming variational Monte Carlo simulations based on 
these projected wavefunctions, we obtain the phase dia- 
gram Fig. [T] and static correlation functions in the spin 
nematic wavefunction. 

We first argue how these mean-field pairing states are 
deformed, when external Zeeman field is applied. A di- 
rect minimization of the mean-field energy dictates that 
all the c?-vectors in these pairing states are restricted 
within a plane perpendicular to the applied field. This 
arrangement makes these pairing states invariant under 
the spin 7r-rotation around the magnetic field. Owing to 
this TT-rotational symmetry, the corresponding projected 
BCS wavefunctions acquire the 'spin-nematic' character; 
ordering of (the transverse component of) the quadrupole 
moments without any spontaneous ordering of (the trans- 
verse) spin moments nor any magnetic crystallization. 
We also show that the 'flat-band' state, which achieves 
the lowest mean-field energy in the strong ferromagnetic 
regime [ Ji [ J2 , actually ends up in the trivial fully po- 
larized state, when projected onto the spin Hilbert space. 

To examine a possible realization of quantum spin ne- 
matic states in the present spin model, we study the en- 
ergetics of the projected (i) Z2 planar state, (ii) Z2 polar 
state, and (iii) SU{2) chiral p-wave state. We focus es- 
pecially on the intermediate coupling regime, where the 
ferromagnetic exchange Ji is about twice as large as the 
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antifcrromagnctic exchange J2, J2 — 0.5|Ji|. Based on 
the variational Monte Carlo analysis, we argued that, in 
a finite range of this intermediate coupling regime, the 
projected Z2 planar state achieves the best optimized 
energy, compared with the energies of other competing 
phases, such as the ferromagnetic phase and coUinear an- 
tifcrromagnctic phase. (See Fig. [T]) We also prove that 
this projected Z2 planar state is accompanied by the 'c?- 
wave' spatial configuration of the quadrupole moments. 
This feature of the wavefunction including the irreducible 
representations under the symmetry group turns out to 
show a perfect agreement with the nature of the quan- 
tum spin nematic phase suggested by Shannon et al^ 
from the exact diagonalization study. 

Motivated by this coincidence, we further calculate the 
spin-spin correlation function of the projected Z2 planar 
state, so as to obtain the static spin structure factor in 
the spin nematic phase. The structure factor thus cal- 
culated exhibits two prominent peaks at the wavevec- 
tors k = (tt, 0) and (0,7r), which signifies the presence 
of strong collinear antiferromagnetic fluctuation. The fi- 
nite size scaling of the peak height concludes that the 
state does not possess any sublattice magnetization in 
the thermodynamic limit, which is consistent with the 
spin-nematic feature of the Z2 planar state. This anti- 
ferromagnetic fluctuation is reminiscent of the neighbor- 
ing collinear antiferromagnetic phase, whereas ferromag- 
netic fluctuation, which is also expected to appear from 
the other neighboring ferromagnetic phase, is completely 
suppressed. 

These observations indicate that the static spin-spin 
correlation function by itself hardly distinguishes the 



current quantum spin nematic phase from the neigh- 
boring collinear antiferromagnetic phase. Indeed, the 
recent exact diagonalization study by Richter et al^ 
reported that the present J1-J2 model exhibits only a 
strong collinear antiferromagnetic correlation in the in- 
termediate coupling regime, J2 ~ 0.5|Ji|. We expect 
that the dynamical magnetic properties in combination 
with these static physical properties could distinguish the 
present spin nematic phase from the collinear antiferro- 
magnetic phase. That is, unlike in the collinear anti- 
ferromagnetic phase, all the gapless Goldstone modes in 
the spin-nematic phase are expected to lose their spectral 
weight in the dynamical spin structure factor in the low- 
energy limit. In fact, a recent calculation based on the 
random phase approximation shows that this is indeed 
the case in the Z2 planar statci^ This feature could be 
sharply contrasted to the dynamical magnetic property 
in the collinear antiferromagnetic phase, where the gap- 
less mode at the (0,7r)-point or (7r,0)-point have a finite 
spectral weight even in the low-energy limit. 
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